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Abstract 

This paper proposes a blind technique for MIMO cognitive radio Secondary Users (SU) to transmit 
in the same band simultaneously with a Primary User (PU) under a maximum interference constraint. In 
the proposed technique, the SU is able to meet the interference constraint of the PU without explicitly 
estimating the interference channel matrix to the PU and without burdening the PU with any interaction 
with the SU. The only condition required of the PU is that for a short time interval it uses a power 
control scheme such that its transmitted power is a monotonic function of the interference inflicted by 
the SU. During this time interval, the SU iteratively modifies the spatial orientation of its transmitted 
signal and measures the effect of this modification on the PU's total transmit power. The entire process 
is based on energy measurements which is very desirable from an implementation point of view. 

I. Introduction 

The emergence of Multiple Input Multiple Output (MIMO) communication opens new di- 
rections and possibilities for Cognitive Radio (CR) networks ^MM- In particular, in underlay 
CR networks, MIMO technology enables the SU to transmit a significant amount of power 
simultaneously in the same band with the PU without interfering with him if the SU utilizes 
separate spatial dimensions than the PU. This spatial separation requires that the interference 
channel from the SU to the PU be known to the SU. Thus, acquiring this knowledge, or operating 
without it, is a major issue of active research [|4|48|]. We consider MIMO primary and secondary 
systems defined as follows: we assume a flat-fading MIMO channel with one primary and 

This work is supported by the ONR under grant N000140910072P00006, the AFOSR under grant FA9550-08- 1-0480, and the 
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multiple SUs. Let Hij G C"^^"^ be the channel matrix between user j's transmitter and the 
PU's receiver. In the underlay CR paradigm, SUs are constrained not to exceed a maximum 
interference level at the PU, i.e. 

\\ii,j^M\'<v, VjVi, (1) 

where is SU j's (j > 1) transmit signal and 77 > is the maximum interference constraint. If 
?7 = 0, the SUs are strictly constrained to transmit only within the null space of the matrix Hy. 

The optimal power allocation for the case of a single SU who knows the matrix H21 in 
addition to its own Channel State Information (CSI) was derived by Zahng and Liang [|ll]. 
For the case of multiple SUs, Scutari at al. yD formulated a competitive game between the 
secondary users. Assuming that the interference matrix to the PU is known by each SU, they 
derived conditions for the existence and uniqueness of a Nash Equilibrium point to the game. 
Zhang et al. [|50 were the first to take into consideration the fact that the interference matrix 
H12 may not be perfectly known (but is partially known) to the SU. They proposed Robust 
Beamforming to assure compliance with the interference constraint of the PU while maximizing 
the SU's throughput. Other work for the case of ^n unknown interference channel with known 
probability distribution is due to Zhang and So [6], who optimized the SU throughput under a 
constraint on the maximum probability that the interference to the PU is above a threshold. 

The underlay concept of CR in general and MIMO CR in particular is that the SU must be 
able to mitigate the interference to the PU blindly without any cooperation. Yi [[SO proposed 
a solution in the case where the SU learns the channel matrix based on channel reciprocity 
between the PU where the SU listens to the PU transmitted signal and estimates Hio's null 
space from the signal's second order statistics. This work was enhanced by Chen et. al. [7]. 
Both works requires channel reciprocity and therefore are restricted to a PU that uses Time 
Division Duplexing (TDD). Once the SU obtains the null space of H12, it does not interfere 
with the PU as long as his signal occupies that null space. 

Other than in the channel reciprocity case, obtaining the value of Hij by the SUs (i.e. the 
interference channel to the PU) requires cooperation from the PU in the estimation phase, e.g. 
where the SU transmits a training sequence, from which the PU estimates Hij and feeds it back 
to the SU. Cooperation of this nature increases system complexity overhead, since it requires a 
handshake between both systems and in addition, the PU needs to be synchronized to the SU's 
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training sequence. This is one of the major technical obstacles that prevents underlay CRs from 
being widespread. 

The objective of this work is to design a simple procedure such that MIMO underlay SU 
can meet the interference constraint to the PU without explicitly estimating the matrix Hi^ and 
without burdening the PU with any handshaking, estimation or synchronization associated with 
the SUs. In the proposed scheme (see Fig. [I]) the PU is not cooperating at all with the SU 
and operates as if it is the only system in the medium (as current PUs operate today). The 
only condition required is that for some short time interval (that may be much shorter than the 
entire learning process), the PU will be using a power control scheme such that its transmitted 
power is a monotonic function of the interference inflicted by the SU. Under this condition, we 
propose a learning algorithm in which the SU is gradually reducing the interference to the PU 
by iteratively modifying the spatial orientation of its transmitted signal and measuring the effect 
of this modification on the PU's total transmit power. The entire process is based on energy 
measurements and on detecting energy variations. Therefore, it does not require any handshake 
or synchronization between the PU and the SU. 

The paper is organized as follows: Section HI] provides the system description and some 
notation. Section |lll] presents a blind approach for realizing the Cyclic Jacobi technique for 
calculating the Eigenvalue Decomposition (EVD) of an unobserved matrix; this will be the 
building block of the blind null space learning algorithm presented in Section |IV] and analyzed 
in Section |Vl Simulations and Conclusions are presented in Sections |VI] and IVIIl respectively. 

IL Problem Formulation 

Consider a flat fading MIMO interference channel with a single PU and a single SU without 
interference cancellation, i.e. each system treats the other system's signal as noise. The PU's 
received signal is given by 

yi (t) = Hnxi (t) + H12X2 (t) + vi (t) , t G N (2) 

where vi(t) is a zero mean stationary noise. In this paper all vectors are column vectors. Let 
A be an / X m complex matrix, then, its null space is defined as A/'(A) = {y e C'" : Ay = 0} 
where = [0, ...,0]^. For simplicity H12 will be denoted by H and the matrix H*H will be 
denoted by G. The time line N will be divided into A^-length intervals referred to as transmission 
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Fig. 1. The addressed cognitive radio scheme. The matrix Hi, 2 is unlcnown to the secondary transmitter and vi (t) is a 
stationary noise (which may include stationary interference). The interference from the SU, Hi2X2(t), is treated as noise, i.e. 
no interference cancellation is performed. The SU measures the energy variations in y2{t) according to which it varies its 
transmission scheme until the interference to the PU becomes sufficiently small such that it does not affect y2{t). 



cycles where, for each cycle, the SU's signal will be constant (this is required only during the 
learning process), i.e. 

M{n-1)N + N') = ^2{{n-1)N+1) 

= ■■■ = x2(A^n + N' -1)= i(n), 

where the time interval nN < t < nN + N' — 1 {N' « N) is the snapshot in which the SU 
measures a function 

Nn+N'-l 

^(^) = iv^ E iiy^wf (4) 

t=Nn 

where y2 is the observed signal at the secondary transmitter that includes the primary system's 
transmitted signal (see Figure [1]). The SU's objective is to learn A/'(H) from {x(n), 
This learning process is carried out in learning stages where each stage consists of K transmission 
cycles. We will index each learning stage by k. The indexing method is depicted in Figure [2l 
During the learning process, the SU varies the interference to the PU by transmitting a different 
interfering signal x(?t,). The secondary transmitter measures y2{t) from which it extracts q{n) 
in order to monitor the PU's transmitted power, i.e. Yl!t=nN ^ Each transmission 

cycle, N, corresponds to the PU's power control cycle, i.e. to the time interval between two 
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-«th transmission cycle- 



[<-sensing snapshot->|<— 



-x,(0 = x(«)- 

zr^ — - 



-►j<-sensing snapshot->| 



t=(n-l)N: 
thePU 
modifies its 
power. 



(= (n-I)N+N': 
the SU 
modifies its 
power. 



time 
unit 



(= iiN: tlie 
PU modifies 
its power. 



t (time) 



;= nN+N': 
tfieSU 
modifies its 
power. 



Fig. 2. The time indexing that is used in this paper, t indexes the basic time unit (pulse time) where N time units constitute 
a transmission cycle that is indexed by ?i. Furthermore, K transmission cycles constitute a learning phase (this is not illustrated 
in this figure). 



consecutive power adaptations made by the PU. In fact the actual q{n) is not important as long 
as it satisfies the following condition. 

Assumption 1: Let k E N, then for every K{k — 1) < n < kK the function q{n) satisfies 
q{n) = /fc(||Hx(ra) IP) where : ]R_|_ — t- R_|_ is a sequence of strictly monotonous continuous 
increasing (decreasing) functions. 

Without loss of generality, we assume that fk is a sequence of monotone increasing functions. 
The most important trait (that will be used in this paper) that follows from Assumption [T] is that 
for every K{k - 1) < n,m < kK, /fe(||Hx(n)|p) > /fe(||Hx(m) p) implies that ||Hx(n)p > 
||Hx(m)|p. The problem is illustrated in Figure [3l From a system point of view, Assumption 
[U means that between two consecutive transmission cycles the primary transmitter's power can 
be modified only due to variations in the interference level from the SU that is at the learning 
process and not in a steady state. 

In the following we provide an example for conditions under which Assumption [T] is satisfied. 
To simplify the exposition, we replace the function g in dH) with 

Nn+N'-l 

^(^) = ]^ E E{||y2(t)f} (5) 

t=nN 

In doing so, we ignore the measurement noise at the secondary transmitter. Such g(n) satisfies 
Assumption [U if for example 
1) The PU's signal has the following form: xi(t) = A/pi(?^)Pis(t) where Pi is a constant 

pre-coding matrix, E{s(t)s*(t)} = I, and pi{n) is the PU's power level at the n}"^ cycle. 

This is satisfied for example if the PU is using beamforming (in this case s(t) is a scalar 
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and Pi is a rank 1 matrix), or using a uniform power allocation (in this case Pi = I) or if 
the PU has a single antenna at the transmitter. 

2) If in addition to condition 1 the PU has an SNR constraint at the receiver or a constant 
rate constraint. Than pi{n) will be a monotonic increasing function of its total noise plus 
interference. 

3) Assume that condition 1 is satisfied and in addition the PU is using OFDMA in which the 
current channel is just one of the bins and assume that the PU is using the water filling 
rule. In that case, pi{n) will be a decreasing function of the interference plus noise. 

Under such conditions and assuming that H21 is unchanged between two transmission cycles. 



The learning process is carried out as follows: At the first cycle (n = 1), the SU transmits 



a low-power signal x(l), such that the interference constraint ([I]) is satisfieqj and measures the 
PU's transmit energy q{l). At the next cycle, the SU transmits X2(2) and measures g(2) and so 
on. Section HVl describes algorithms for the SU to learn the null space of the interference channel 
matrix H based on these measurements, i.e. to approximate A/'(H) from {x(?7,), q{n)}^^^ where 
the accuracy can be arbitrary small for sufficiently large T. The algorithm is not limited to 
networks with a single SU; it is also valid for networks with multiple SUs as long as only one 
system modifies its power allocation during its learning process. This fact enables a new SU to 
join the channel in which multiple SUs coexist with the PU in a steady-state (i.e. a case where 
each SU meets its interference constraint given in ©). 

III. Blind Jacobi Diagonalization 

In this section we present a blind approach for realizing the Jacobi technique for calculating the 
EVD of an unobserved Hermitian matrix G assuming that only S{G,x) = x*Gx is observed. 

' This can be obtained by transmitting a very low-power signal at the first cycle. If qp (n) is not affected, the power of x can 
be gradually increased until q{n) is affected. 



we have. 



q,{n) = E{\\y,{t)r} 

= pi(r2)E{H2iPis(t)s(t)*H2iPl} + E{||v2(t)f } 
= pi(n)Tr{H2iPiH2iPl} + E{||v2(t)f } 
= aipi(n) + 02 



(6) 
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This algorithm will be the building block of the blind null space learning algorithm that is 
presented in Section |IVl 

A. The Jacobi Technique 

The Jacobi technique obtains the EVD of the Hermitian matrix G via a series of 2-dimensional 
rotation that eliminates two off-diagonal elements in each phase (indexed by k). It begins with 
setting Aq = G and then performing the following rotation operations 

Afc+i = V,AfcV*,fc = l,2... (7) 

where = R;^m(6', 0) is an Ut x Ut unitary rotation matrix whose p, q entry is given by: 

cos(6') if p = q E {m, 1} 

g-«0 gin(^) if p = m ^ q = I 
[Rz,m(^, = \ -e'^ sin(^) if p = l^q = m (8) 

1 if p = q ^ {m, 1} 

otherwise 
For each k, the vales of 6, cp are chosen such that [Afc+i]^ ,„ = 0, in words, 6 and </> are chosen to 
zero Afe's /,m and m, / off diagonal entries. The values of /, m are chosen in step k according 
to a function J : N — > x {l,...,nt} i.e Jk = {Ik.'mk). It is the choice of Jk that 

differs between different Jacobi techniques. For example, in the classic Jacobi technique, the off 
diagonal elements are chosen according to 

(lk,mk) = argmax (|[Afc]z,m|) (9) 

{{l,m.):l>m} 

which corresponds the maximal off-diagonal entrjo. In the cyclic Jacobi method the rotation rule 
is defined as follows: 

Definition 1: Jk = {lk,rnk) is a function such that 1 < Ik < Ut — 1 and Ik < rrik < rit where 
each pair (/,m) is chosen once in each cycle. Unless otherwise stated it is assumed that Ik < Iq 
if k < q and that rrik < niq if k < q and Ik = Iq. 

For example if Ut = 3 then Ji = (1, 2), J2 = (1, 3), J3 = (2, 3), J4 = (1, 2).... 

^ Recall that A is Hermitian therefore it is sufficient to restrict m > I. 
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x(n) e C- 



H 



q{n) e 



/.(•) 



Fig. 3. Block Diagram of the Blind Null Space Learning Problem. The SU's objective is to learn the null space of H by inserting 
a series of {x(n)}„gN and measuring q{n) as output. The only information that can be extracted is that ||Hx(n)||^ > ||Hx(/)|p 
if q{n) > q{l) for every (fc - 1)K <l,n< kK where fc G N. 



B. Blind Jacobi Step 

In the proposed cognitive radio scenario, if the SU wishes to perform the Jacobi step it has 
to do so without observing G and but observing only /fc(||Hx(n) |p) = /fc(x*(?7,)Gx(?7,)) as 
depicted in Figure [H The following proposition shows how to do this: 

Theorem 1: Let G be an rit x rit Hermitian matrix, let m(^) 0) be defined in Eq. ([8]), and 
let ^(G, a) = a*Ga. The values of 6 and (p that eliminate the /, m entry of G i.e. 

[R,,^(^,0)GR,V(^,0)],,^ = O (10) 

are given by 

(^,0) =argexte,^^(G,r,,„(^,0)) (11) 

where argext denotes an extreme point and ri^(6',0) denotes the /th column of R;,„(6',0). 
Furthermore, every local minimum point of the function 5 (G, r/ ^(^5 0)) is also an absolute 
minimum point., i.e. , let T* = {(^*, 0*) = 7* G : 3 e > , ^ (G r;,„,(7*)) < S (G, r;,„,(7)) , 
V II7-7II < e}' then 5(G, ri,^(7i)) = S{G,vi^rn{l2))^ '^li,l2 ^ T*. The same statement 
applies to local maxima. 

Proof: See Appendix El 
Theorem [T] asserts that the Jacobi step can be carried out by a blind two-dimensional opti- 
mization in which every local minimum is a global minimum. This is a very important property 
since it is easy to identify if an optimization algorithm has converged to a local minimum. Note 
that Theorem [H applies also if ^(Gjx) is not observed but /^(^(Gjx) is observed instead. 
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C. Reducing the Two-Dimensional Optimization into Two One -Dimensional Optimizations 

Although the Jacobi step can be implemented blindly, it requires a two dimensional optimiza- 
tion over the parameters Q and 0. This may be very difficult in practice since each search point is 
obtained via a transmission cycle. Fortunately, this two dimensional optimization can be carried 
out optimality by two one-dimensional optimizations as shown in following theorem: 

Theorem 2: Let S'(G, x) = x*Gx and let ri.m{9, 0) be Ri,m(^, 0)'s /"^ column where R/,m(^, 4>] 
is defined in ([8]). The optimal Jacobi parameters 9 and in (fTTI) can be achieved by 

= argminS'(G,rK7r/3,0)) (12) 

(Ae[-7r,7r] 

{ e if - 7r/4 < ^ < 7r/4 

I 9 — sign(6')7r/2 otherwise 

where sign(6') = 1 if 6* > and —1 otherwise and 

e~ = arg min ^ (c, r,,„(e, 0)) (14) 

Proof: See Appendix lAl 
Comment: Note that the rotation angle 9^ is restricted to the interval [— 7r/4, — '/r/4]. In Section 
lIV-Bl it is shown that this restriction guarantees a globally linear convergence rate and ultimately 
a quadratic convergence rate. 



be carried out using line searches. This is 



In practice, the optimizations in (fT2l) . (fT4l) will 
because, the function S'(G,x) will not be observec^but only /fc(S'(G,x)) will be observed (see 
Figure [3l). The complexity of the one dimensional line search in (fT2l) can be drastically reduced if 
one is looking for a minimum (or maximum) and if the objective function has a single minimum 
(local and global). Under these conditions a binary search can be invoked. This is important since 
it reduces the complexity of the line search drastically from 0(l/?7) to 0(log(l/?7)) where rj is 
the line search accuracy. In the sequel, it will be shown how to solve both (fT2l) and (fT4l) using 
a binary search. We begin with the optimization in (fT4)) which can be written as 

0fc = arg min /fc (5 (G, rz,„(7r/3, 0))) (15) 

0e[-7r,7r] 

'in it is siiown that 5(G,x) is known, the problem can be simplified drastically where G can be obtained precisely by 
a finite number of transmission cycles. 
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Note that 

S{G,Ti^^rn{0,4>)) = 008^(9) \gi^i\ + sin^(6') \gm,m\ 

(16) 

- \gi^rn\ sm{2e) cos(0 + Zgi^rn) 
It follows that during each learning phase, that is for every (A; — 1)K < n < kK, the function 
fk iS{G, ri,^(0.5, 0))) is equivalent to 

Wk{(p) = fk{A + B cos{(f) + Zgi^J) (17) 

where /fc is a monotone function and A > \B\. 

Proposition 3: Consider the function w in (flTI) then 

a. is monotonic on [7r/2,7r] and non-monotonic on (0,7r/2) if and only if 

\w{0) -w{7r/2)\<\w{7r/2) -w{7r)\ (18) 

b. Let G (— vr, tt) be a minimumQ point of w{(f)) then: 

i. Assume ([l8]) is true, then, e [0,7r/2] if w{tv/2) > w{tt) and G [-tt, -7r/2] if 
w{tx /2) > w{n). 

ii. Assume ([l8]) is false, then, G [-7r/2,0] if w{'k/2) > w{0) and G [vr/2,7r] if 
M;(7r/2) < w{0). 

Proof: This follows immediately from the fact that A cos(0 — + i? is affine symmetriclfl 
in the intervals „i — 7r/2, Zgim + 7r/2] and [Z(7i,.m + 7r/2, /lgi,m + 37r/2] with unique extreme 
point. □ 
Proposition [3] can determine, via 3 points, a 7r/2-length closed-interval in which = Zf?/,™ is 
the unique local (and therefore global) minimum. It is now possible to invoke a binary search to 
approximate Zgi^rn by where for any accuracy, say 77 > 0, it takes an additional — log2(2?7/7r) 
points (transmission cycles) to ensure G [Zgi^rn—f]^ ^gi,m+v] - exactly the same manner it can 
be shown that the optimal value of 9 in (fT2)) can be approximated 77 closely using 3 — log2(?7/7r) 
transmission cycles. The algorithm is summarized in Table HI 

''since w is assumed monotone and continuous and its argument is a 2-k periodical sinusoid sucli point always exist if 

■w{0) 7^ w{n). 

'a function / : R — > R is affine symmetric if there exist some a G R such that f{x — a) = f{a — x). 
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TABLE I 

Line Search 



function: z — LineSearch(w, Zmin, Zmax, ??) 

1) Initialize: z = 2;niax/2; a, &, c = 0; L = Zmax- 

2) If \w (Zmax/2) - 10(0)1 > \w (2:niax/2) - W (^niax)], Set a= I. 

3) If to(0) < li; (2max/ 2), Set 6 = L 

4) Ifw( ) > w (zniax/ 2) set c = L 

5) Zmax = 2maxa(l - &) + 2max(l - a)c; 

7) Wliile (I 

2m ax -2^111 in 

\>v) 

8) 2 = ( + 

9) If (^(2111111) < iu(2max)), Set Zmax = Z, Otfierwise, set Zmin = z. 

10) end while 
end LineSearch 



IV. Blind Null Space learning via the cyclic Jacobi Technique 

In this section we present a Blind Null Space Learning (BNSL) algorithm for solving the 
blind null space learning problem based on Theorems [U- 12] and on the line search algorithm in 
Table HI The SU is using a pre-coding matrix that is being updated at each learning stage, let 
Wfc be that pre-coding matrix at stage k. 

A. The Blind Null Space Learning (BNSL) Algorithm 

Let USV* be H's Singular Value Decomposition, where V and U are rit x rit and rir x Ur 
unitary matrices respectively and assume that nt> rir. The matrix S is an rirXrit diagonal matrix 
with real nonnegative diagonal entries ai, that are arranged as ai > cr2, > ■ ■ ■ > 0"^ > 0. 
We assume with no loss of generality that = d{= Rank(H)). Note that 

A/'(H) = span (v„^+i, v„J (19) 

where Vj denotes V's ith column. From the SU point of view, it is sufficient to learn A/'(G) 
(which is equal to A/'(H) because G = VAV* where A = S^S). It is possible to diagonalize 
G via the blind Jacobi algorithm given in Table HIl The secondary user's initial (i.e. at A; = 0) 
pre-coding matrix is Wq = I. Each Jacobi step is equivalent to a learning stage (indexed by 
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k, see Figure [21) and is composed of K transmission cycles where the secondary user obtains 
{9k,(l)k) at its end, and updates its pre-coding matrix 

Wk = Wk-iRi„m,{ekAk) (20) 

Obtaining {9k, 4>k) at each learning stage requires two binary searches where the value of fk{S{G, 
WfeF;^ 0(n)))) for each search point (6'„,0„) is obtained by a transmission cycle in 

which the SU is transmitting 

X2(t) = i(n) = Wfc„ir,,,„,(^„,0„) e C"* ^^^^ 

, V(n -1)N <t <nN 

i.e. for each k, the SU performs two one-dimensional binary searches to obtain {6k, (pk) (the 
minimum of fk{S{G, Wk — ii"ife,mfe (^5 0)))) from the set \^ici^^m^{6n, <pn ), fk{S{G, Wfe„ir/^ 

,mk i^n 

, 4'n)))}n={k-i)K- algorithm proceeds, 6k approaches zero (this will be shown on Section 

IVT) . a fact that can be used for a stopping criteria. 

Assume that the BNSL algorithm is performed until k = kg. Then the SU pre-coding matrix 
Tk, is given by 



Tfc = [wf%...,wf= 1 (22) 



where wf is W^'s ith column and ii,i2, ■■■,int is an indexing such that 

(wJO*Gw|^ < (w£0*Gw|^ < ■ ■ ■ < (w^;)*Gwt^ (23) 
Thus, the interference power that the SU inflict to the PU is bounded 

||Hi2Xi2f <P2||Hi2TfcJ|2 (24) 

where p2 is the SU's transmit power. In Section |V] it will be shown that ||Hi2TfcP converges 
quadratically to zero as is bounded by for sufficiently small 7] and ultimately bounded by 0{Tf). 

It is important to note the eigenvalues of G cannot be obtained by the BNSL algorithm. This 
leaves the SU with the problem of how to determine which of the columns of W^,.^ corresponds 
to Hi2's null space, if it doesn't know its's rank in advance. The problem however can be solved 
due to the fact that v, G N{G) iff S{G, v^) = S{G, 2vi). Denote hk : C"' M+ as 

hk{-^) = fk{S{G,^)), (25) 

thus, by setting x(?i) = Vj and x(n + 1) = 2vj it follows that Vj G A/'(G) iff hk{St{n)) = 
hk{x{n + 1)). The same approximately applies for W^:^ when kg is sufficiently large. 
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TABLE 11 

The Blind Null Space Learning (BNSL) Algorithm 



Given the sequence {/ifcjfegN defined in i25l 
function: W = BNSL({/ifc}fegN, n*) 

1) Initialization: k = 1, Wo = I„t, Aj = 2^,\/j < 0. 

2) wiiile (maxjg{fe_„(„_i)/2,...,fc} > ^). 

3) w{(f>) = hk (Wk-iri^^jn^{Tv/3,(j})). 

4) set 0fe = LineSearch(M, 0, n, rj) 

5) to(6() = /ifc (Wfc_lr,,,,„J6(,,^fe))• 
6) set Ok = LineSearch(TO, 0, 7r/2, 77) 

7) set 9k according to (T3j. 

8) set Afc = l^fel 

9) set Wfc = Wfc_iR,,,„,(4,0fc), W = Wfc 

10) k^k + l. 

11) end while 



5. The Reduced Complexity Blind Null Space Learning (RC-BNSL) Algorithm 

In the BNSL algorithm, each step requires two line searches where each search point is 
obtained by a transmission cycle. These transmission cycles are the dominant latency factor in 
the learning process since the rest of the calculations are performed off line at the secondary 
device processing unit. Roughly speaking, the complexity of the cyclic Jacobi technique grows 
like nf which is the dimension of the matrix G (the convergence and complexity will be discussed 
in Section |Vl). In this Section we present an algorithm that converts the problem of blind null 
space learning of the Ut x ut matrix G into an equivalent problem of blind null Space learning of 
Ut — rir matrices where each is an x riy. matrix. This is possible due to the fact that G = H*H 
is a rif-rank matrix and therefore has a — dimensional null spaced. The resulting complexity 
grows like {rit — nr)n^, therefore, the algorithm is efficient if rit is sufficiently larger than Ur 
(which is a practical case since SU systems will be more sophisticated than the PU). 
The idea behind the RC-BNSL algorithm is described in the following observation: 
Observation 4: Let H e C"'^^"* be an n^-rank (nt > rir) matrix and let U = [ui, u„,,+i] G 

*If the matrix G was known, the fact that rank(G) = Ur could have been utilized by QR decomposition (which cannot be 
done blindly) prior to the diagonalization. 
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Qntxinr+1) orthonormal matrix (that is, a matrix whose columns are an orthonormal set 

i.e. U*U = I) and let H = HU. If Ci G AfiU) then u = Ufl G Af{U). 
Proof: This is due to Hu = HUu = since u G A/'(HU). 
The RC-BNSL algorithm is carried out as follows: The secondary user begins with an initial 
pre-coding matrix U'^°) G C"r>i{nt~nr) ^\Yich is composed of the last Ut — riy. columns of some 
unitary matrix W G C"*''"'. Let H^^ = HU(°) G C"'''("'+^), then there exists at least one 



degree of freedom in this channel. The SU can apply the BNSL algorithm on G^^) = H^q^'H 



r(i) 



and obtains a pre-coding matrix \J^^^ such that U*^^) = linifc^oo U^^^ (in Section |V] it is shown 



(1) 



that the limit exists) and that 







A 







(1) 











(26) 



Now that G*^^^ is diagonalized the first degree of freedom is given by v*^^^ = U'-'^^uf ^ where 
u^^-* is the first column of \J^^\ (that lies in the null space of G*^^^). The SU then can gain 
an additional degree of freedom by applying the BNSL algorithm on the following [rii + 1) 
equivalent channel 

Hg) = HU(^) (27) 

where U*^^-' G C'^'^"'^"'"^ is obtained by concatenating the — 1 column of the initial unitary 
matrix W with the last column of U^^^ multiplied by U^°\ i.e. let tj'^^^ = [u2^\ u^^^^^^^] 
then U*^^) = [wr„^_i, U'^^^U'^^^]. This equivalent channel is then diagonalized using the BNSL 
algorithm to obtain U'^^\ We now have two degrees of freedom given by v^^^ = U'-^^uf ''and v'^^\ 
This process is repeated until all W column are used. The RC-BNSL algorithm is summarized 
in Table Unl 



V. Convergence 

In this section we study the convergence and complexity properties of the BNSL algorithm. 
Recall that this algorithm is in fact a blind implementation of the cyclic Jacobi technique whose 
convergence properties have been extensively studied over the last 50 years. However, while the 
convergence properties of the Cyclic Jacobi technique directly apply to its Blind implementation 
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TABLE III 

The Reduced Complexity Blind Null Space Learning (RC-BNSL) Algorithm 



Let W G C"'^"' be a unitary matrix, 
function: [vi ■ • • v„j_„^] = RC_BNSL({/ifc}fcgN, w^, ^t) 
Initialize U*"^ = [wnj , w„j]. 
for m = 1, ... ,nt — rir 

a(x) = /ifc(U(™-i)x) 
U(m) ^ BNSL(a,n^ + 1) 

u™ = [fir\...,Qi:ii] 

u(™' = [wi7!„_„,u(™-i)u('")] 

end for 



due to Theorem [T] they cannot be applied directly to the BNSL algorithm. This is due to fact 
that in the latter algorithm the rotation angles 9^, (pk are approximated via a line search of finite 
accuracy rj for every k (see Table |Ill lines ll® while in the previous, 9^, (pk, are equivalent 
(according to Theorem ^ to the rotation angles of the Cyclic Jacobi technique. Moreover, we 
would like to make this line search accuracy as small as possible (that is, to make t] as large 
possible) in order to reduce the number of transmission cycles. It is therefore very important to 
understand how r] affects the performance of the BNSL algorithm. In this section we will extend 
the classic convergence results of the Cyclic Jacobi technique to the BNSL algorithm and show 
what is the required accuracy in the line search that assures convergence and a minimal level of 
interference inflicted by the SU to the PU. 

A. Overview of the previous work on the convergence of the Jacobi technique 

The convergence of the Jacobi technique has been studied extensively over the last sixty years. 
The first convergence proof of the Cyclic Jacobi technique for complex Hermitian matrices was 
given by Foster and Henrici |1(|], which proved that if in each step k: 1) The off diagonal entry 
satisfies laf^^^fcl < ^'^I'^t^mfel' where af'^^ is A^'s /,m entry (defined in (|7])) and < Cjt < 6 < 1 
where h is independent of k. 2) The rotation angle (in this paper it is 9^ defined in (O) lies in 
some close interval A E (— 7r/2, 7r/2), then the cyclic technique converges to G's EVD. These 
conditions are satisfied by the BNSL algorithm for any line search accuracy for which 9k ^ 0. 
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Thus, to guarantee convergence of the BNSL algorithm to G's EVD, one should follow the 
following rule: If 9^ = 0, phase k should be repeated with r] being decreased so that 9^ ^ 0. 
If any further decrease in t] does not change 9^ ^ 0, then af^ is already zero. This however 
does not indicate of the convergence rate and does not take into considerations that cannot be 
infinitely decreased. 



Henrici, and Zimmermann llllh proved that the cyclic Jacobi algorithm for real symmetric 
matrices has a global linear convergence rateO that depends on the matrix size rit if the rotation 
angle 9^ G [— 7r/4, 7r/4] for every k. Fernando [12] extended this result to complex Hermitian 
matrices. A very important result is the ultimate quadratic convergencetl rate of the Cyclic Jacobi 
technique that was shown by Henrici [[ibI for complex Hermitian matrices with well separated 
eigenvalues and later enhanced by Wilkinson lll4ll . The most recent and comprehensive result 
of the quadratic convergence of Cyclic Jacobi technique that includes multiple and clusters 
of eigenvalues is due to Hari fl5^. Once the Cyclic Jacobi algorithm reaches its quadratic 
convergence rate it takes a very small number of cycles to reach any desirable accuracy, however, 
there is no rigorous bound on the number of the required cycles to reach that rate. Brent and 
Luk [16] have argued heuristically that this number is 0{log{nt)) cycles for nt x rit matrices. 
This seems to be the case in practice 11171 page 429] where such a rapid decline is obtained after 
three to four cycles [18:, page 197]. For further reading the reader is referred to lll7l-ll9n. 



B. Global Linear Convergence rate 



The global linear convergence of the Cyclic Jacobi technique was derived by Fernando il2|] : 
Theorem 5 (nlW.Theorem 4): Let G be a finite dimensional nxn complex Hermitian matrix 
and Pk denote the norm of the off diagonal upper triangular (or lower triangular) part of (as 
defined in (|7])) and let m = n(n — l)/2 (which is a cycle length). Then 

P.+„.<pP.,P=(l-2-(--i)(--^)/^)'/' 



if 



TT 



max {%\} < - 

v<k+m 4 



(28) 



(29) 



'a sequence a„ is said to have a linear convergence rate of < /3 < 1 if ja„+i| < /3ja„|. 

*A sequence is said to have a quadratic convergence rate if there exist some /3 > such that |an+i| < /3|anp. 
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Although Theorem \5\ is the tightest closed form bound on the global convergence rate of the 
cyclic Jacobi technique, the practical convergence rate is much faster as discussed in Section 

The condition of Theorem \5\ is satisfied in the BNSL algorithm as we show in Theorem [2l 
However, as described in Table |IIl the angles 9^ and (p^ are approximated via a line search. 
Thus the off diagonal elements are not completely annihilated, i.e. [Afc];^. ~ instead of 
[Ak]i^,mk = 0. In the following Theorem we show what is the required line search accuracy that 
guarantees the convergence rate in (|28l) . 

Theorem 6: Let G be a finite dimensional rit x rit complex Hermitian matrix and Pk denote 
the norm of the off diagonal upper triangular (or lower triangular) part of = W^_^GWfc_i 
where is defined in (|20|) (see also Table|Ill) and let m = nt{nt — l)/2. Let rj be the accuracy 
of the line search, then the BNSL algorithm has a globally linear convergence rate that is given 

< Pi (1 - 2^("*-2)(-'-i)/2) + [nf - nt){7 + 2V2W\\Gf (30) 

Proof: See Appendix |B] 
To demonstrate this result we substitute rj = aPgm/\\G\\ and obtain 

1 /2 

Pk+m <Pk{l- 2^(2-nOK-i) ^ ^ 2v^) (n^ - n^)) (31) 

It follows that for the BNSL to have a linear convergence rate, it is sufficient that the accuracy 
be at least: 



n 


3 


4 


5 


6 


8 


77 X 


8 X 10-2 


2 X 10^2 


7 X 10-3 


1 X 10-3 


2 X 10-^ 



(32) 



Note that the linear convergence coefficient in Theorem [6] may be very small for a large 
number of transmitting antennas Ut. This may be very bad if this bound turns out to be tight, 
i.e. if p is very close to the true convergence coefficient. In the sequel it will be shown that the 
actual convergence rate of the BNSL is much faster than the bound in Theorem [6l 

C. An Asymptotic Quadratic Convergence Rate 

So far it has been shown that for a right chose of 77, the BNSL algorithm converges and that 
for sufficiently small rj it has at least a global linear convergence rate. In the following theorem 
and corollary it is shown what is the desired rj to guarantee an asymptotic quadratic convergence 
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rate and what 77 guarantees that the SU will meet the maximal interference constraint of the PU. 

Theorem 7: Let rj be the accuracy of the line search, {A^}"^^ be G's eigenvalues and let 

35 = min |Ai — Aj.| (33) 

Let Pk be the norm of the off diagonal upper triangular part of = W^_^GWfc_i where 

is defined in (l20l) (see also Table |Ill) and let m = nt{nt — l)/2. Assume that the BNSL algorithm 

has reached a stage where P| < 5^/8, then 

+2 {2ntnr — nl — nr) ?7^||G|p 

Proof: See Appendix O 
Theorem |7] shows that to guarantee quadratic convergence rate the accuracy should be much 
smaller than P^, that is, let /cq be an integer such that P^^ < 5^/8 then 

ft„,„<o((^)'j (35) 

if ?7 << P|y. This implies that once P^ becomes very small such that P^ « 7], one cannot 
guarantee that Pk+m will be smaller than P| but only be smaller than 0(77). A fact that motivates 
derivation of a bound on the interference power that the SU inflict to the PU as a function of 

V- 

Corollary 8: Let be the SU's pre-coding matrix defined in (l22l) . Assume that the conditions 
of Theorem |7] are satisfied, then 

m^^Tk+mf < O j + O (^^^ j + 2 {2ntnr - nl - n,) TylGf/^ (36) 

Proof: This is an immediate consequence of Theorem [TT] in Appendix O 
Note that the quantity ||Hi2Tfc+m|P is the only one that interest the SU that applies the 
BNSL algorithm. Furthermore, once r] >> P^ the dominant term in the (l36l) will be 0{r]) i.e., 
the interference power to the PU will approximately satisfy 

\\B.uTm+kf < 2 {2ntnr - nl - n,) v^\\G\\yS (37) 
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This allows choosing 77 to guarantee a maximum interference level to the PU, an observation 
that will be very useful in the simulation part. Theorem |7] and Corollary [8] also imply that the 
line search accuracy need not be constant during the entire BNSL algorithm but can be refined 
as the algorithm goes on. 

The asymptomatic quadratic convergence rate of Theorem |7] and Corollary [8] is determined by 
1/5 where 35 is the minimal gap between G's eigenvalues. In addition, the quadratic convergence 
rate takes effect only after P| < 5/8. Such a condition implies that if 5 is very small, it will take 
the BNSL many cycles to reach its quadratic convergence rate. This is problematic since MIMO 
wireless channels may have very close singular values (recall that H12 square singular values are 
equal to G's first rir eigenvalues). If we were using the optimal Cyclic Jacobi techniq ue ( i.e. no 
errors due to finite line search accuracy) this would not have a practical implications lllSn since 
quadratic decrease in Pk that is independent of 5 occurs prior to the phase where P| < 5/8. In 
the following theorem we extend this result to the BNSL algorithm. 

Theorem 9: Let 77 be the accuracy of the line search, {Xi}]^^^ be G's eigenvalues such that there 
exist a cluster of eigenvalues, i.e. \i = \ + ^i, for I = nt — v + 1, rit where YMlnt-v+i ^ ^■ 
Assume that the rest of the non equal eigenvalues are well separated, i.e. 6c >> where 

35c = min(Ai U A2) (38) 

Ai = {|Ai -Xr\:l<l<r<nt-v, Xi ^ A,.} ^^^^ 
A2 = {\Xi-X\:l<l<nt-v} 

Then, once the BNSL algorithm reaches a stage such that 25c ^^^^ < Pk < ^c/^^ and if 
Tj « Pi, then 

P.+.<0((^)') (40) 

Proof: See Appendix |Dj 
In the presence of very close eigenvectors cluster, i.e. << ^c, the distance 5c will 

be much greater than 5. In that case, a quadratic decrease in will occur even before P^ 
becomes smaller than 6/^/8 but only satisfies 25ca/^7^ ^ Pk ^ This quadratic decrease 
brakes down and become slower as P| becomes smaller than 25c a/X]/ ih From a practical point 
of view, this is not a problem if one is not interested in decreasing P| more than 25c ^Xl; ^1 
which may be very small. Nevertheless, Pk will eventually decrease quadratically as P| becomes 
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Fig. 4. Simulation results for different line search accuracy values 77 of the BNSL algorithm for obtaining the null space of H 
with Tit = 3 transmitting antennas and Ur = 2 antennas at the PU receiver. The vertical axis represents the square of the sum 
of the magnitude of the upper off-diagonal entries of G = H*H while the horizontal axis represent the number of complete 
cycles of the BNSL algorithm, i.e. [ni — nt)/2 learning phases. The matrix G where normalized such that ||G||^ = 1. We used 
200 Monte-Carlo trails where the entries of H are i.i.d. complex Gaussian Random variables. 



smaller than 5/8 as required by Theorem |7l This phenomena is for the Cyclic Jacobi technique 



in yj 



VL SIMULATIONS 

Figure |4] presents simulation results of the BNSL algorithm for different levels of line search 
accuracies. It is shown that for sufficiently small rj the algorithm converges quadratically. The 
quadratic decrease breaks down when the value of becomes as small as an order of magnitude 
of rj. This result is consistent with the bound that is proposed in Theorem |7l In addition, Figure 
|4] shows that as long as r] is smaller or equal to the decrease in Pk+m is almost not affected 
by different line search accuracies as demonstrated by the fact that the decrease is approximately 
the same for 77 = — 10, —20, —40 at the first cycle as well as for 77 = —20, —40, at the second 
cycle. The simulation shows that this phenomena, which is consistent with the asymptomatic 
behavior as indicated by Theorem |7l is also valid before the algorithm reaches its asymptomatic 
behavior. Figure [5] describe the interference decrease as a function of the transmission cycles 
for different line search accuracies. The result is consistent with Corollary [8] where the ultimate 
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Fig. 5. Simulation results for adoptive line search accuracy values of the BNSL algorithm for obtaining the null space of H 
with nt — 3 transmitting antennas and different number of antennas at the PU receiver rir. The vertical axis represents a bound 
on the norm of the interference to the PU while the horizontal axis represent the number of transmission cycles. The matrix G 
where normalized such that ||G||^ — 1. We used 100 Monte-Carlo trails where the entries of H are i.i.d. complex Gaussian 
Random variables. 



level of the interference is bounded by 0{rf). Similar to Figured the interference decrease in 
Figure [5] is approximately the same for rj = —10,-20,-40 at the first cycle as well as for 
7] = —20, —40, at the second cycle. Moreover, Figure [5] shows that increasing the values of rj 
reduces the number of transmission cycles drastically as discussed in Section IIII-CI 

Both Figures |4] and \5\ suggest that using a low line search accuracy (i.e. larger 77) in the first 
cycle and increasing it from one cycle to the other may reduce the overall transmission cycles 
with no significant performance loss. This idea is put into practice in Figure |6] where we present 
simulation results of the interference decrease to the PU as a function of the transmission cycles 
for an increasing line search accuracy. 

VII. CONCLUSIONS 

We proposed a blind technique for interference mitigation by secondary cognitive users based 
on a blind implementation of the well known cyclic Jacobi technique. The only condition required 
is that the primary user will be using a power control scheme such that for a short time interval, 
its transmitted power be effected monotonically by the interference from one SU. This includes 
also a case where there are multiple secondary users in a steady state, i.e. their interference power 
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Fig. 6. Simulation results of BNSL algorithm for non constant line search accuracy and different numbers of transmitting 
antenna at the SU transmitter. The line search accuracy ri is -6dB, -8dB -15dB at the first, second and third cycles respectively. 
The PU has rir — 2 transmitting antennas. The vertical axis represents the reduction in the interference to the PU while the 
horizontal axis represent the number of transmission cycles. We used 200 Monte-Carlo trails where the entries of H are i.i.d. 
complex Gaussian Random variables. 



is constant during this time period. The entire learning process is based on energy measurements 
and on detecting energy variations. This means that the secondary users are not required to be 
synchronized to the PU pulse time. Furthermore the proposed learning scheme is independent 
of the PU transmission scheme (i.e. coding, modulation) as long as the power allocation is a 
monotonous function of the interference from the secondary user. 

The Convergence properties of the BNSL algorithm were also explored in this paper. It was 
shown that the BNSL algorithm maintains the convergence properties of the Cyclic Jacobi 
technique, that is, a global linear and an asymptotically quadratic convergence rates, as long 
as the line search accuracy is sufficiently small. It was also shown in simulation that just like in 
the cyclic Jacobi technique, the BNSL algorithm reaches its quadratic convergence rate in just 
three to four cycles. Furthermore, we obtained a bound on the interference that the SU inflicts 
to the PU as a function of the line search accuracy of the BNSL algorithm and provided a 
mechanism for choosing this line search accuracy to reduce the number of transmission cycles 
while maintaining low levels of interference to the PU. 

It is important to stress that the BNSL algorithm is not necessarily limited to energy measure- 
ments taken by the SU and to PU that apply power control. The BNSL algorithm can also be 
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implemented in other scenarios as long as the SU can learn whether the interference it inflicts to 
the PU has increased or decreased between one transmission cycle to the other. For example, the 
secondary user can learn about the interference by observing the PU's modulation scheme since 
this is too a function of the the interference at the PU. An alternative way for implementing 
the BNSL algorithm is for the SU to listen to the PU's control channel and extract information 
about the interference it inflict to the PU. 



Appendix A 
Proof of Theorems [Hand [2] 
The idea behind the proof is that each /, m rotation is equivalent to diagonalizing a 2 x 2 
matrix Gi^m using the rotation matrix R(^, 0) as follows 



gi,i gi,m 

9m,l 9m,m 



R(^, 



-i(f> 



sm{9) 



-e'^sm{9) cos(^) 



(41) 



Let ri{6, 0), r2{9, 0) be R(6', 0)'s first and second columns respectively. It follows that 



)GRr 



)Gr,. 



)G; mfi(6', I 



The first part of the Theorem [H follows directly from the Rayleigh-Ritz Theorem [see e.g. 
Theorem 4.2.2] that asserts that 



(42) 



20L 



.A 



l,m 



min x*G/,rn.x 



where i? = {x G : ||x|| = 1} and Xf^^ is G/^m's minimal eigenvalue. It follows that 



\ mm 



)G;mri( 



where 



arg min f^(^, < 



)G^mfl(^,' 



(43) 



(44) 



(45) 



ax 



Equality ([441) is due to the fact that for every x G -B there exist 6*, G M such that f 1(6', 
where,|a| = l,a G C and because for every G G C^^^ and x G the function S satisfies 
S'(G,ax) = ^(Gjx) it follows that max^gg 5(6/,^, x) < max^^^gK ^(G/^rn., ? 0))- In addi- 
tion, because f 1(6', (j)) E B for every 9, (j) we have that max^.^^ S{Gi^rn-, x) > maxg^^gK S{Gi^rn-, 
f 1(6', (j))) which establishes (|44l) . Note that fi(6', 0) is the eigenvector that corresponds A™™ and 
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since ri and r2 are orthonormal, it follows that r2 is the eigenvector that corresponds to Gi,m's 
maximal eigenvalue Xf^^. Hence, 



R(^,0)G,^R*( 



\ mm fi 



A 



max 



(46) 



and because [Rz,m(^, 0)GR;,„(^, 0)];,^ = [R(^, 0)Gi,^R*(^, 0)]i,2 = 0, the desired result 
follows. 

It remains to show the second part of the theorem. The objective function for minimization is 
S (G, v^e, 0)) = [R,,„(^, <p)GRUe, = cos^e) \g^\ ^^^^ 

9m,m\ \9i,m\ sin(2^) cos(0 + /-gi Yn) 
which is a continuous and differentiable (of any order) function of 6 and 0. Recall that G > 0, 
thus gii are real non-negative numbers and therefore will be written as gu instead of \gii\. We 
first assume that gi2 7^ 0. Setting the gradient to zero yields 



= -2^1,2 cos(2^) cos {^gi^m + 0) + gi,i{- sin(2e)) + ^2,2 sin(2^) 
= 5(1,2 sin(26') sin (Z^j^^ + 0) 
and the solution is 



(48) 



hn \ ( Tib 71 

-l)'^^o + Y,avr-Z(7i,„j , _, -712 + - + Tra j (49) 



(50) 



where 

Itan-if-^) if gu^g„.„. 

f if 9ll = 9mm 

We begin with the family of suspected points ((— l)"6'o + ^, an — Zgi^rn)- Since ^(G, ri{9, 0)) = 
S{G, ri{6 + TT, (p + 27r)), \/9, G M it is sufficient to investigate the following subset 

{9, 0) e {7i = (^0, -Z5/,m), 72 = (^0 + 7r/2, -Zgi^J 

, 73 = (-^0, vr - Zfi'/.m), 74(^/2 - 6*0, vr - Zgi^rn} 
Furthermore, because ri(7]^) = — ri(73) ^"1(72) = ~i"i(74) sufficient to check the 
points 7i and 72. To investigate these points we calculate the Hassian of ^(G, ri(9, cp)) that is 
given by 

y^{S)efi = 4 sin(26') \gi^rn\ cos (0 + Zgi^rn) + 2 cos(26') {gm.^ra - 9l,l) 

V^{S)e,4> = 2 cos(2e) |^/,^| sin (0 + Zgi^^) (52) 

V^('S')<^,<^ = sin(26') \gi^ra\ cos (0 + Zfi-/,™) 
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and its primary minor (determinant) 

det {S7\S)) = \gi^m\ (2 \gi,m\ (cos (20 + 2Z^i,„) - cos(4^)) 
+ sin(46') {gm,m - 9i,i) cos (0 + Zfif/.m)) 
is equal to 4|(7imp for 7j,j = 1,...,4. Therefore these are local extreme points. In order to 
determine which is the maximum and which is the minimum we will substitute these points in 

|cos(eo)|2(4gf 2+(gi,i-g2,2)^) , 

45fi2 otherwise 



2(gi,i-g2,2) j£ / 

V'{S)e,oiOo + ^/2,-Zg,,m) = { 1™^'°' yn ^ y22 ^^^^ 

—4(712 otherwise 



Note that (1541) and (1551) are of opposite signs, therefore, one of which corresponds to a local 
minimum while the other to a local maximum. It remains to check the family of suspected 
points (y, — 7i2 + f + ^'^)- By substituting it in (|53] ) we see that it is equal to -4|(7; for every 
a,b E 7j. It follows that these suspected points are not extreme points. This establishes that for 
fi'1,2 7^ all of the of local minimum points (which are infinitely countable) of S{G,ri(9, (p)) 
are equivalent, i.e. for every local minimum points {9,(p), {0' such that (^,0) ^ {d',<P') we 
have ^(G, ri(^, 0)) = S'(G, ri(6'', 0')). Thus, every local minimum is also a global minimum. 
In the case where gi^2 = 0, the target function is 

5(G,r;,„(^,0)) = cos2(^) \g^\+sm\e) \g^^^\ (56) 

which obviously satisfies the conditions of Theorem [TJ This establishes the proof of Theorem [T] 
To proof Theorem |2] we substitute 6 = n/S in (l47l) 



S (G, r;,m(7r/3, 0)) = ^3 \gi^rn\ cos (0 + Zgi^jri) + — + ^^"^'""^ (57) 

It follows that in the blind one-dimensional minimization min,^ S* (G, r/ m(7r/3, 0)), cos(0 + 
'^9i,m) = 1 which is obtained by = —Zgi m since the line search is carried out in the 
interval G [— 7r,7r]. By performing blindly the optimization in (fT4l) we obtain the minimum 
that corresponds to either 7^ or 72 in (|5TI) . Then, the value of 9 is chosen such that 9 lies in 
the interval [— 7r/4, 7r/4], thus (6', 0) G {71,73} where 71,73 are defined in (|5TI) . □ 
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Appendix B 
Proof of Theorem [6] 

Consider the first cycle of the cyclic Jacobi technique, i.e. k = 1,2, ...,nt{nt — l)/2. Denote 
the number of rotated elements in the Ith row hy bi = rit — I and let 

E 1=1 = {2n, - 1 - 1)1/2 



Cl 



k\l,j+l\ 



w{i, k) = z{j, k) 



(58) 



Note that W{Q, k) = P|. In every cycle, each entry is eliminated once, we therefore denote A^'s 
p, q entry before its annihilation as gq^p{t) where t denotes the number of changes since A; = 0. 
After gq,p{t) is annihilated, it will be denoted by gq,p{i) where i is the number of changes after 
the annihilation. The diagonal entries of A^ will be denoted by x since we are not interested in 
their values in the course of the proof. This is illustrated in the following example of a 4 x 4 
matrix 



Ao = G 

/ ^?i,i(0) 91,2(0) (71,3(0) (71,4(0) \ 

^72,l(0) (72,2 ( 0) (72,3 ( 0) (72,4 ( 0) 

^73,l(0) (73,2 ( 0) (73,3 ( 0) (73,4 ( 0) 

V ^?4,l(0) (74,2(0) (74,3(0) (74,4(0) / 
A2 

^ X e (71,4(2) \ 

^2,i(0) X 5-2,3(2) 52,4(1) 

e 53,2(2) X 53,4(1) 

V 54,1(2) 54,2(1) 54,3(1) X ) 



Ai 



X e 
e X 
53,1(1) 53,2(1) 



51,3(1) 5m(1) \ 
52,3(1) 52,4(1) 
X 53,4(0) 



X 



\ 54,1(1) 54,2(1) 54,3(0) 
A3 

/ X 51,2(1) 51,3(0) e \ 

52,1(1) X 52,3(2) 52,4(2 

53,1(0) 53,2(2) X 53,4(2 

e 54,2(2) 54,3(2) X 



(59) 



For arbitrary rit, after the first ci sweeps Ac^'s first column is equal to the following vector: 

[x, h^iijit - 3), 5„,_i,i(0), e]"^ (60) 

and 



^(l,ci) < |52,i(rit-3)p + ... + |5n,-i,i(0)P + |e| 



(61) 
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For q = 2, rit we have 

ggAnt - q - I) < cos(6'„,_i)^g,i(ra-g-2) -e^'^"t-i5(g,„,(l)sin(6'„,„i) 

(62) 

9,,iil) < cos(Vi)^.,i(0)-e*'^^(7,,g+2(l)sm(e3) 
^,,i(0) < |e|cos(^,)-e*S,,,+i(l)sin(^g) 

where 1) = e. Bounds on {^g,i(0}r=o'^^ '^^'^ obtained recursively (i.e. obtaining bound 

on i(0), substituting and obtaining bound on and so on...) 

~g,,^in, -q-l)<e UZ,' cos (9^) - EJl"' e^"^^ sin (9,) g,,^^{l) UT'li {9.) ^^^^ 
= v(g)Tz(g) + |6|nr;cos(e.) 

where v, y e C"'-" such that [v(g)], = -e^*^+i(7,,,+,(l) and [y(g)],- = sin (^,+,_i) n"l^+, cos (^„)- 
It follows that 

\~g,,,{n, - g - l)p < |y^(g)v(g)p + |6p ^os^ (9.,) 



<iiy(9)iniv(g)ip + iepnr;cos2(^„) 



(64) 



Proposition 10: 

\\y{q)r = l-l[cos{9,) (65) 

i=<j 

Proof: This is shown by induction. Note that 

n— 1 rt— 1 

\\y{q)r = Y,^m\9,) J] cos^l^.) (66) 

i=g v=i+l 

where 

m 

JJq = l,if 1 > m. (67) 

i=i 

Assume that 

m—l m—1 m—1 

1 - J] cos'(^,) = sin'(^i) n ^°^'(^-) (6^) 

1=5 1=5 t>=i+l 



is true for some m G N. Then, for m + 1 we have 



+1 

YZ~" sm^(^0 nr=.+i cos2(^.) + sin2(^„) Y[Z^^, cos\9.) (69) 

m—1 • 2rn \ TT™^1 



cos2(^m) sin^(^.) nrii ^os^l^.) + sin^ 
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where the last equality is due ( |67] ). According to the supposition (|68] ) 

= i-nr=,cos2(^.) 

By substituting Proposition \T0\ into (|641 ) one obtains 

i=l / v=q 
' V ' ^ ^ ' 



(70) 

□ 



(71) 



= Z{q,Ci_) <1 

thus, 

|^,,i(ni-g-l)|2< M - J] cos2(^,) J Z(g,ci) + |e|2 (72) 
and by summing both sides of (|72l) over g = 2, 

^(i,ci) < e:i2 (i - nri,cos2(^o) + (n, - i)i6p 



9=2 ^ 

< (1 - rcico+2Cos2(^.))W^(0,0) + (n, - l)|6p 



(73) 



where the last inequality is due to Pc.= Ci) + Z(l, ci), Vr(0, 0) = Pq, and because Pk is 
a monotonically decreasing sequence It follows that 

Z(l, ci) = sin^ (^,o+2,cJ W^(0, 0) + (rii - l)|ep (74) 

where 

Sin^ (^e,_,+2,cj = 1 - n ^°''(^~^) (^5) 

j=Ci_i+2 

and |6'i| < 16*4 1. Thus, 

Pei = ci) + ci) < 1^(0, 0) = Po (76) 



substituting (1741) we obtain 



ci) < 1^(0, 0) cos^ (^2,cJ - {nt - l)|ep (77) 



'Forsythe and Henrici llOll showed that the sequence Pk is a monotonically decreasing sequence. 
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Now that this relation is established, it will be applied to Ac/s lower {rit — 1) x {rit — 1) 
block-diagonal. To do that we use the fact that for Jci+i = (2, 3), = (2,^^ — 1) the 
sum of squares of the first column (which is equal to the first row) remains unchanged, i.e. 
Yl^Li |[Afc]i,„P is constant for k = Ci,Ci + 1,...,C2. Thus 

W{l,ci) < W{1 - cos^ (*c,_,+2,c,) - {rit - IM (78) 

Continuing this way we obtain 

I I I 

W{1, q) < W{0, 0) n cos^ (vl>,^_,+2,c,) -e'Y, U cos' (vl>,„_,+2,c.) (79) 

j=l j=l v=j+l 

thus 

Z{l,ci) = sin^ (^e,_,+2,cj W{1 - l,Q„i) + {nt - l)\e'\ 

< W{0, 0) sin^ (*,_,H-2,.) Uy=\ cos^ (*c,_,+2,c,) (80) 

After a complete cycle 

< W^(0, 0) Eri7' (*e,_,+2,c,) n;=i cos^ (^c,-,H-2,c,) (81) 
Similar to proposition [TOl it can be shown that 

n l~\ n 

J2 sin' (Ti) n cos^ ir,) = ^11 ^os' (r,) (82) 
1=1 j=i j=i 

Thus 

PL. < WiO, 0) (l - n;=1 COS^ (^c,_.H-2,c,)) 

From (1751) we have 

cki 

COS^ (^c,_,+2,c,) > n ^°^'(^-) (^"^^ 
?;=£,_ 1+2 

and therefore 

PL. < w{o, 0) (i - n;if ntc,_,H-2 cos^i^.)) 



(83) 



(85) 
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Recall that \9i\ < n/i, therefore 



Pl_, < Vr(0,0) (l-2-("-2){-i)/2) 



[Ei=i Ei=i(^ - j)2V-'"+^+9«-45 - j (86) 

< 1^(0) (1 - 2-("-2){n-l)/2) ^ |e2| 

It remains to relate e to the accuracy of the line search -q. Note that the error e in ( f86l) results 
from the two finite- accuracy (of rj accuracy) line-searches in Table |Il] on lines |4] and [6l If r/ were 
zero, Afc's /, m off diagonal entry would be zero after the A;th sweep, i.e. 

M(rP*,0°P*) = (87) 

where 



u{e,<P) = |[R,,^(e,0)AfcR* (e,0)],,^|2 = 4{af)' sin' (7,,™ + 



(88) 



+ (2cos(2^)af^„cos (Zaf^„ + 0) + sin(2^) (af,, - a^J) 

and iO'k'^ ■, (j^k^) is the value given in Theorem [2] when substituting G = A^. Let {6k,(f)k) be 
the non optimal value that is obtain by the two line searches (on Line |4] and Line (6]), then 

|ep = maxfc n(6'fc, 0fc). The error u{9k,(j)k) can be bounded because 0^^* = Zaf„, thus 0^ = 
-Zo;^^ + ?70 where \r]^\ < 77. Thus 

Uiiek,<Pk) = 4(af,j2sin2 (7,,^ + 0,) < 4(af^)V < 2||G||r/2 (89) 

The second term 

U2iek, <Pk) = {2al^ cos {ri^) 008^(2^^) + sin(2^fc) {a^ - al^^)f (90) 

Note that if = a^^, the value 9k = 0°^^ G 0,7r/4 since the line search will not miss these 
points. Now for the case where af; 7^ aj^^ we have 9k = 91 + rje where 

91 = ]- tan-' (xk) (91) 



2 



and 



2|af,rr^|cos(r/<^) 



Note that 

U2{9k, 0fc) = (2 cos (r/^) af;„ (cos (291) - 2r/, sin (2^*)) 
+ Kz - <™) (sin (2^^) + 2 cos (2^) ve))' 



(93) 
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where {9*, (jf) is a point on the line that connects the points (6'^^*, 4>°k^), {Ok, <Pk)- By substituting 
(1911) we obtain 

MOM = ( '-M<-^t'-'''=<- - 4ry, sin (2^*) cos (v^) a^^ + 2ry, cos (26*) {a^, - a^J 

(94) 



by (1921 ) and by the fact that sinusoidal is bounded by one and by \rig\ < r] we obtain 

U2{eM < W (2| sin(2^*)|a^^ + cos(2^*) |a^, - a^jf 

< V (4sin2(2r)|af,^p+2sin(4e*)|af,^||4 - a^J + cos2(2e*)|4 - alJ") (95) 

< V (2|af,^p + 2sin(4r)|af;„J|< - a^J + |< - a^J^ + 2|af;„j2) 

u2{ekAk) < W (2||Gf + v^||G||||G|| + \\Gf) (96) 

Thus 

|ep = maxfeu(4,</)fc) <2(7 + 2v^)r/2||G||2 (97) 
This expression is substituted into (|86l) and the desired result follows. 

Appendix C 
Proof of Theorem [7] 

We first assume that G's eigenvalues are all distinct. From (l63l) it follows that 
Similarly to the derivation of (|72|) we have 

(99) 



and by summing both sides of (|99l) (similarly to the derivation of (|73l) ) over g = 2, ...,nf it 
follows that 




(100) 

< (E;i2'sin^(^.)) mO,0) + K-l)|ep 
Now that we established this relation we can apply it to the reduced Ut — I + \ lower block 
diagonal and obtain 

Z{l,ci) < (E;U_,+iSm^(%)) W^(0) + K-OkP (101) 
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After a complete cycle we have 



, (102) 

< w{o, 0) e;^'"'^/' sin2(%) + k'l Eri7'(^* - 

We now relate ^"l^f "^^/%m2(^^) to W{0,0). Let Ai,...,A„, be G's eigenvalues and let 

35 = min |A; — Ami (103) 
and assume that the algorithm has reached a phase where 

p2 = y^/^o, k) < 6^8 (104) 

,fc „fc 12 |„fc \ „fe _|_\ _|_ \ \ |2\| \ \ 12 \^k \|2 \ 12 

Hi '^mml ~ "mm "T -r /\; Am | ^ I^V; Am\ l^ii \^mm j 



furthermore, [ | 13u Theorem 1] we have |afi - Ai| < 5/2. Thus, 

|af,-a^„| >25-5/2-5/2 = 5 (105) 

Recall that the optimal rotation angle satisfies tan(20^'") = 2|af^^^ |/|af^;^ — am^mfcl while the 
actual the rotation angel is 9k = 0'^^ + rjg. It follows that 

I sin2(4)| < I sm\el'' )| + \r^gsm{2eT)\ < i|2^,P 
+ \rie\tan{29T)\ < ^tan\2eT ) + \ve\ tan(2^f )| (106) 



< g2 + 2|?7g| ''g < ''ga" H 



It follows that 



(107) 



^m(m-l)/2 . 2(n ^ ^ ^m(m-l)/2 / ^.mf 2M^W(0,k)\ 

= 3^W^(O,A:) + ^^^Mz!^0y(O;i) 
By substituting (|107l) into (11021) one obtains 

< W^(0,0) [^W{0,0) + + ^ K -n,) , (108) 

It remains to relate 77^ to the accuracy of the line search r]. As a result of the error on LineH 
r]g depends on 77^ as well. Form the proof of Theorem \T\ we know that if an accurate line search 
were invoked, it would produce (pk = —^<^im.- However, due to the finite accuracy 77, the line 
search yields 0^ = —Zaf^ + r]^, where |?7(^| < r/. Thus, 9^ is obtained by optimizing a slightly 
perturbed version of w{9), say w{9), due to the substitution of (pk into hk (See Table HI] Line [5]) 
i.e. 

w{9) = S{Ak, ri,mi9, = hk{cos\9)al - cos{r]^) sm{29)al^ + sm''{9)ai^J (109) 
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We first assume that af; 7^ CLmm- If both line searches were accurate, the optimal value of 9 
would be 



9T = \te.n-'{pk) (110) 



2\a'' I 

where pk = — — ^3-^- If one takes into consideration the non-optimality of the line-search on 



Line |4] and ignores the non-optimality of the line search on Line [6] then w{9) is given in (11091) 
and the optimal value would be 



91 = ^tan ^ {pkCos{r]^)) (111) 



The difference \9l^' - 9l\ is 



l^r* - ^fcl = 1 1 tan ^ {pk cos(?7^)) - ^ tan ^ (p^ 



^ sin(>?;)pfe| ^ 2 bfcl 



(112) 



+1 



where It]^] < i]^. It can be easily that 



< 1 ^ (113) 



cos^{r]^)pl + l |cos(r/^)| 
Because 6'^. = 91 + rifp and l?]^! < 77, the accumulated effect of the finite accuracy on both Line 
in and Line [6] is bounded by 

2 

Ve<V + r^^ (114) 
I cos(?7j| 

Assuming that 77 is sufficiently small, (e.g. i] < n/20) we obtain 

Ve<Qv/5 (115) 
By substituting (11151) and ([97]) into (11081) it follows that 



(116) 



Pi,., < w{o, 0) (;iiy(o, 0) + v^P^Vwm] 

+ (10 + 2v^)(n?-n0r/2||G|P 

Thus, as long as 77 is smaller than W{0, 0), the BNSL will have a quadratic convergence rate for 
G that does not have multiple eigenvalues, i.e. all eigenvalues are distinct. This is not sufficient 
since we are interested in a matrix G rit — rir with zero eigenvalues. 

To extend the proof to the case where the matrix G has Ut — rir zero eigenvalues and 
distinct eigenvalues we shall use the following theorem: 
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A 



(117) 



Theorem 11 (hi 8(1 Theorem 9.5.1): Let A be an rit x rij Hermitian matrix with eigenvalues 
{AJJ!!^ that satisfy Ai 7^ A2 7^ ■ ■ ■ 7^ 7^ A„^+i = A„,,+2 = ■ ■ ■ = A^^ = A. Consider the 
following partition: 

' Ai B 
B A2 

where Ai is ra,. x ra,. and A2 is (ut — rij) x [rit — rir) and let 5' > 0. If ||(Ai — AI)"^|| < 1/5', 
then 

IIA2-AIII < ||B||V(5' (118) 

We now apply Theorem [TTI to the BNSL Algorithm. Let A5^,A2,B'^ be A^'s submatrices 
that correspond to the partition in (II 171) . Recall that in our case, A = 0, thus, (11051) implies that 
||A^|| > 6. Furthermore, by nlUi Corollary 6.3.4], the matrix A^' is invertible, thus ||(A^')~^|| < 
1/6, and form Theorem [TTI it follows that 



lA^II < WBkf/S 



(119) 



To show how this theorem leads to quadratic convergence we first show that once the BNSL 
algorithm reaches a stage where (11041) is satisfied, the affiliation of the diagonal entries in the 
upper Ur X rir -block remains unchanged, i.e. if A^ satisfies (11041) then 



argmin \Xi — afi\ = argmin |Am — aj^|, VL C {1, n^} 



(120) 



Note that 



< sin^ {Ok) (2cos(^fc)a;';^^^ cos (0^ - Zaf^^^J + sin {Ok) - «^„„^J)' 

(121) 

and that for every 9k such that 1 < A; < c„^ (11061) is satisfied. Thus 



< sin^ m (at™, + sin (^,) (af^,^. - a^,^,^J) ^ 

< sin^ {Ok) (at., + [{almS/^' + ^Ve al„ Js) 

< sin^Ok) {al^^ + {Sy 4 + 26voS/2)y')' 

< sinM^.)(at„,^+5(l/4 + r/,)V2)^ 

< f (1 + Ar]e) (1/2 + ^1/2 + 7^0 



(122) 
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Note that for 77 < 1/100 and considering (11151) it follows that 

Ka 0-655 (123) 

which establishes (11201) . Now that (11201) is established, it follows that for every l,m such that 

I 7^ m and 1 < / < n^, we have \afi — a^^\ > 6. 
After Cnt^r rotations (11021) can be written as 

< EriT' z{i, c„j + 16^| Er=iK - ^^24) 

< W{0, 0) + Eri;;+i ^(^, c„J + le^l j:Z,{n, - I) 

Recall that |ep < maxku{6k, (pk) where Ui{9k,(t)k) and U2{9k,4>k) are defined in ( f89l ) and (l90l) . 
From (|95l) we have 

M2(^fc,0fc) < V(4P,' + 4Pfc||G|| + ||G||2) (125) 
Because |af; — a^^l > 5, (11071) is satisfied and similarly to (II 161) we obtain 

< 1^(0,0) (^^l^(O,O) + r/^^0l^(O:i)) 
+2 {2ntnr - - n,) r/2(9iy(o, 0)/2 + 4^1^(0, 0)||G|| + ||Gf) (126) 

It remains to bound the term Yll^d +1 ^(^' ^nr)- Note that for every 6k such that 1 < A; < c„^, 
( 11061 ) is satisfied. Let Q = {{l,m) : 1 < I < n,. < m < Ut] . Because af^ and are located in 
A^^ and Ag respectively, it follows that for every k such that {lk,mk) E Q 

|ag+\,|' < K^f + sm\9k)\alJ^ for < g < 



(127) 



and from (11061) 



la'= P la* I 

I .mi. I , f-v I li. .mi. I 



< K,/ + ( + 2r/e^^ ) 14,/ for mfc < g < 

la 



I' < |a,^,„J' + ( %^ + 2r/,%^ ) lafj', for < g < 



(128) 



These can be bounded by 

l„fc+i 12 i^fe+i 12 

Thus, for every k < 



<W\0,0){l + jj) + ^W^/\0,0) (129) 



ErC.^(^,c.J= "e E |ay^<0((™^l)')+0((^^^^ (130) 

q=n^+l t=q+l \^ ^ J ^ ^ 

February 3, 2012 DRAFT 



36 



This, together with (11261 ) and ( I115D show 

2" 



^n — r — 



-0 



W(0,0) 



r)2tyi/2(0,0) 



o 



TiW'-^/^ (0,0) 



2 {2ntnr 



nz 



rir) rf\\G\ 



Since Pk is a decreasing sequence, the desired result follows. 



(131) 



Appendix D 
Proof of Theorem [9] 

Let VfcAV^ = Afc be A^'s EVD, and let 

A'^ = VfcAV^ 



where 



A^ = VfcAV^ 
A = diag(Ai, ■ ■ ■ , A„i_^_r,^j^^,^Xj^^) 

r V 

A = diag(Oj^,ei,--- 

nt—v 



(132) 



(133) 



Similarly to the derivation of (11051) . there exists a permutation such that |a{} — a^^.\ > 6c for 
all /, m. such that / < and l<l<nt — v — rorrit — v — r<l<nt — v + 1 and I > rit — v. 
For such a permutation, let A;,, be partitioned as 



A A: \k \k 

^11 "^12 -^13 

\k \k \k 

^21 ^22 ^23 

\k \k \k 

-^31 -^32 -^33 



where A^2 ^ C^^' and A^g e Then by [UJ, Lemma 2.3] we have that 



(134) 



p2 

K\\oS< :r^,for / = 2,3. 



26. 



(135) 



where ||AfJ|ofr represent the sum of squares of Af/s off diagonal entries. The rest of the proof 
is identical the proof of Theorem |7] from equation (11191) and forward since (|119|) is satisfied by 
setting 



\k _ Ak -ofe _ r A fc \k 



13J' "^2 



A 2 \k 

^22 ^23 

Afc Afc 

^32 "^33 



(136) 
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